ctis_smoothed  %>%
  filter(keep) %>%
  mutate(state=toupper(state)) %>%
  select(date,
         state, 
         imputed_beta,
         beta_estimate_smoothed, 
         beta_estimate_spline_smoothed) %>%
  pivot_longer(contains("beta")) %>%
  mutate(name = case_when(
    name == "beta_estimate_smoothed" ~ "LOESS smoothed",
    name == "beta_estimate_spline_smoothed" ~ "Spline smoothed",
    name == "imputed_beta" ~ "Survey Value"
  )) %>%
  ggplot(aes(x=date, y=value, color = name, 
             alpha = name, linewidth=name)) +
  geom_line() +
  facet_wrap(~state, ncol=4) +
  scale_alpha_manual(values=c("LOESS smoothed" = .9, 
                            "Spline smoothed" = .9, 
                              "Survey Value"=.3),
                     name='') +
  scale_linewidth_manual(values = c("LOESS smoothed" = 1.05, 
                            "Spline smoothed" = 1.05, 
                              "Survey Value"=.5)) +
  scale_color_manual(values=c("#3381FF", "#B58746", "#26900F"),
                     name ='') +
  guides(alpha="none",
         linewidth="none",
         color = guide_legend(override.aes = list(linewidth = 3,
                                                  alpha = c("LOESS smoothed" = .9, 
                            "Spline smoothed" = .9, 
                              "Survey Value"=.3)),
                            nrow=3)) +
  theme_c(legend.position="top") +
  ylim(0,1) +
  scale_x_date(date_breaks="3 months",
               date_labels = "%b %Y") + 
  labs(y = TeX("Survey Estimate of $\\beta$"),
       x= "",
       title = TeX("Comparing Approaches for Smoothing Survey Estimates of $\\beta$"))

ctis_smoothed  %>%
  mutate(state=toupper(state)) %>%
  select(date,
         state, 
         contains("s_untested")) %>%
  pivot_longer(contains("s_untested")) %>%
  mutate(name = case_when(
    name == "s_untested_smoothed" ~ "LOESS smoothed",
    name == "imputed_s_untested" ~ "Survey Value"
  )) %>%
  ggplot(aes(x=date, y=value, color = name, 
             alpha = name, linewidth=name)) +
  geom_line() +
  facet_wrap(~state, ncol=4) +
  scale_alpha_manual(values=c("LOESS smoothed" = .9, 
                              "Survey Value"=.6),
                     name='') +
  scale_linewidth_manual(values = c("LOESS smoothed" = 1.02, 
                              "Survey Value"=.9)) +
  scale_color_manual(values=c("Survey Value" ="black", "LOESS smoothed"="darkred"),
                     name ='') +
  guides(alpha="none",
         linewidth="none",
         color = guide_legend(override.aes = list(
           linewidth = 3,
           alpha = c("LOESS smoothed" = .9, 
                     "Survey Value"=.3)),
                            nrow=3)) +
  theme_c(legend.position="top") +
  scale_x_date(date_breaks="3 months",
               date_labels = "%b %Y") + 
  labs(y = TeX("Survey Estimate of $Pr(S_1|untested)$"),
       x= "",
       title = TeX("Comparing Approaches for Smoothing Survey Estimates of $Pr(S_1|untested)$"))

# raw ctis data

ctis_raw <- readRDS(here('data/data_raw/ctis_all_states.RDS'))

states_keep <- ctis_smoothed  %>%
  filter(keep) %>% pull(state) %>% unique()

ctis <- ctis_raw %>%
  select(signal, state=geo_value, date, value) %>%
  pivot_wider(names_from=signal, values_from=value) %>%
  filter(state %in% states_keep) %>%
  mutate(beta=smoothed_wscreening_tested_positive_14d/
           smoothed_wtested_positive_14d,
         state=toupper(state))

ctis %>%
  group_by(state) %>%
  summarize(median = median(beta, na.rm=TRUE),
            mean = mean(beta, na.rm=TRUE),
            q1 = quantile(beta, .025,  na.rm=TRUE),
            q2 = quantile(beta, .975,  na.rm=TRUE)) %>%
  ggplot(aes(x= fct_reorder(state, mean, .desc=TRUE),
             y = mean)) +
  geom_errorbar(aes(ymin=q1,ymax=q2),width=.2) +
  geom_point(color="darkred", size=2) +
  scale_y_continuous(n.breaks=10) +
  labs(y= TeX("Empirical Estimates of $\\beta$"),
       x="State",
       title = TeX("Empirical Estimates of $\\beta$ from the COVID-19 Trends and Impact Survey"),
       subtitle = "2.7% Percentile, 97.5% Percentile, and Median By State" ) +
  theme_c(axis.title=element_text(size=12))

ggsave(here("figures/emp-estimate-beta-dist-by-state.pdf"))

#
# Logistic transformation of the Beta CDF.
#
f.beta <- function(alpha, beta, x, lower=0, upper=1) {
  p <- pbeta(x, alpha, beta)
 log(p/(1-p))
}
#
# Sums of squares.
#
delta <- function(fit, actual) sum((fit-actual)^2)
#
# The objective function handles the transformed parameters `theta` and
# uses `f.beta` and `delta` to fit the values and measure their discrepancies.
#
objective <- function(theta, x, prob, ...) {
  ab <- exp(theta) # Parameters are the *logs* of alpha and beta
  fit <- f.beta(ab[1], ab[2], x, ...)
  return (delta(fit, prob))
}
#
# Solve two problems.
#
par(mfrow=c(1,2))
alpha <- 15; beta <- 22 # The true parameters
for (x in list(c(1e-3, 2e-3), c(1/3, 2/3))) {
  x.p <- f.beta(alpha, beta, x)        # The correct values of the p_i
  start <- log(c(1e1, 1e1))            # A good guess is useful here
  sol <- nlm(objective, start, x=x, prob=x.p, lower=0, upper=1,
             typsize=c(1,1), fscale=1e-12, gradtol=1e-12)
  parms <- exp(sol$estimate)           # Estimates of alpha and beta
  #
  # Display the actual and estimated values.
  #
  print(rbind(Actual=c(alpha=alpha, beta=beta), Fit=parms))
  #
  # Plot the true and estimated CDFs.
  #      
  curve(pbeta(x, alpha, beta), 0, 1, n=1001, lwd=2)
  curve(pbeta(x, parms[1], parms[2]), n=1001, add=TRUE, col="Red")
  points(x, pbeta(x, alpha, beta))
}

ctis %>%
  group_by(state) %>%
  summarize(median = median(beta, na.rm=TRUE),
            mean = mean(beta, na.rm=TRUE),
            q1 = quantile(beta, .025,  na.rm=TRUE),
            q2 = quantile(beta, .975,  na.rm=TRUE)) %>%
  ggplot(aes(y=fct_reorder(state, mean, .desc=TRUE))) +
  ggridges::geom_density_ridges(
    aes(x=beta,
        y=fct_reorder(state, beta, .fun =median)),
    data=ctis,
    fill="#596281",
    quantile_lines=TRUE,
    quantiles = c(.025,.975),
    quantile_fun=function(beta,...) c(
                                    #  quantile(beta,.025, na.rm=TRUE),
                                      median(beta,na.rm=TRUE)
                                    #  quantile(beta,.975, na.rm=TRUE)
                                      )) +
  # geom_errorbar(aes( xmin=q1,xmax=q2),width=.2) +
  scale_x_continuous(n.breaks=10, limits=c(0,1)) +
  labs(x= TeX("Empirical Estimates of $\\beta$"),
       y="State",
       title = TeX("Empirical Estimates of $\\beta$ from the COVID-19 Trends and Impact Survey"),
       subtitle = "Density Estimate of Distribution and Median by State" ) +
  theme_c(axis.title=element_text(size=12),
          axis.text.x=element_text(size=11)) 

ctis %>%
  ggplot(aes(x = beta, fill=state)) +
  geom_density(alpha=.6) +
  scale_fill_viridis(option="mako", discrete=TRUE) +
  theme_c()

ctis %>%
  group_by(state) %>%
  summarize(mean=mean(beta,na.rm=TRUE)) %>%
  pull(mean) %>%
  mean()

ctis %>%
  group_by(state) %>%
  summarize(mean=mean(smoothed_wcli,na.rm=TRUE)) %>%
  pull(mean) %>%
  mean()

source(here("R/05-state-analysis.R"))


beta_shape = get_shape(ctis_smoothed,
                             option="beta",
                       quiet=TRUE)




tibble(x = seq(0, 1, length = 10^(4)),
       original = with(prior_params, beta_density(x,
                                     mean = beta_mean,
                                     sd = beta_sd,
                                     bounds = beta_bounds)),
       ctis= dbeta(x, shape1=beta_shape[1], 
                   shape2=beta_shape[2])) %>%
  pivot_longer(c(original,ctis))%>%
   ggplot(aes(x=x, y = value, fill=name)) +
  geom_line(alpha = .8) + 
  geom_ribbon(aes(x=x,ymin=0, ymax=value,fill=name),
              alpha=.6) +
  theme_c(legend.text=element_text(size = 10)) +
  viridis::scale_fill_viridis(discrete=TRUE,
                             option = "rocket", begin=.3,end=.8) +
  labs(x = "Value",
       y = "Probability Density", 
       title = latex2exp::TeX("Definition of Prior for $\\beta$"),
       fill ='')

s_untested_shape = get_shape(ctis_smoothed,
                             option="s_untested",
                       quiet=TRUE)




tibble(x = seq(0, 1, length = 10^(4)),
       original = with(prior_params, beta_density(x,
                                     mean = s_untested_mean,
                                     sd = s_untested_sd,
                                     bounds = s_untested_bounds)),
       ctis= dbeta(x, shape1=s_untested_shape[1], 
                   shape2=s_untested_shape[2])) %>%
  pivot_longer(c(original,ctis))%>%
   ggplot(aes(x=x, y = value, fill=name)) +
  geom_line(alpha = .8) + 
  geom_ribbon(aes(x=x,ymin=0, ymax=value,fill=name),
              alpha=.6) +
  theme_c(legend.text=element_text(size = 10)) +
  viridis::scale_fill_viridis(discrete=TRUE,
                             option = "rocket", begin=.3,end=.8) +
  labs(x = "Value",
       y = "Probability Density", 
       title = latex2exp::TeX("Definition of Prior for $Pr(S_1|untested)$"),
       fill ='')